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In order to reliably compute the longitudinal structure functions in decaying and forced turbu- 
lence, local isotropy is examined with the aid of the isotropic expression of the incompressible con- 
ditions for the second and third order structure functions. Furthermore, the Karman-Howarth- 
Kolmogorov relation is investigated to examine the effects of external forcing and temporally 
decreasing of the second order structure function. On the basis of these investigations, the scal- 
ing range and exponents of the longitudinal structure functions are determined for decaying 
and forced turbulence with the aid of the extended-self-similarity (ESS) method. We find that 
(„'s are smaller, for n > 4, in decaying turbulence than in forced turbulence. The reasons for 
this discrepancy are discussed. Analysis of the local slopes of the structure functions is used to 
justify the ESS method. 
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§1. Introduction 

Study of structure function is useful in the analysis of 
the scaling property of fully-developed turbulence IIP In 
particular, the structure function of the longitudinal ve- 
locity increment .clearly shows how the scaling deviates 
from K41 r scalingcl' as the order of the structure function 
increases.™ The structure function of a transverse veloc- 
ity increment also exhibits deviation from K41 scaling, 
but the values of the scaling exponents appear to dif- 
fer from those of the longitudinal variety, reflcxtin^thc 
influence of vortical structures in turbulence.Q&BB&tf 

All measurements in direct numerical simulations 
(DNS's) and experiments include certain limitations or 
restrictions. Examples of these restrictions are the lim- 
itation of the Reynolds number to finite values, the 
size of the samples, and the degree of deviation from 
isotropy, homogeneity, and stationarity that are allowed 
to the system. Therefore the results obtained from such 
data must be studied carefully and thoughtfully. To 
be more specific, consider the DNS of forced homoge- 
neous isotropic turbulence with a finite Reynolds num- 
ber, in which the random forces are assumed to be sta- 
tistically homogeneous, isotropic and Gaussian with a 
given spectrum. In order to realize a high Reynolds tur- 
bulence number, the force spectrum is assumed to have 
compact support at the very low wavenumber range, 
1 < k < 2 ~ 3. This restriction introduces a problem of 
convergence to large scale values with respect to isotropy, 
homogeneity and stationarity. Since the small number of 
Fourier modes in the forced band have a large amount 
of kinetic energy, and their characteristic times are long, 
fluctuation of these modes significantly affects the statis- 
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tics of large-scale turbulence, which possibly pertain to 
the inertial region. The difficulty lies in determining for 
how long statistical quantities must be averaged to assure 
statistical homogeneity, isotropy and stationarity of the 
system, and in discovering the largest scale not affected 
by the external forces. 

In simulation of decaying turbulence, similar problems 
arise, because the fluctuations of the Fourier modes in 
the lower band are large at t = 0. Hence the issues to 
be addressed are the effect of non-stationarity on the in- 
tensity of fluctuations, the anisotropic effect due to large 
eddies, and the size of the largest scale that is not influ- 
enced by a decrease in the intensity of fluctuations. 

In laboratory tests, such as experiments using the tur- 
bulence in_ a cylinder with disks rotating in opposite 
directionsjHP we find a situation similar to the case of 
forced DNS. The geometry of the experimental appara- 
tus imposes certain limitations on the isotropy and ho- 
mogeneity of the system, although the sample size is rel- 
atively larger than in the DNS case and the Reynolds 
number is relatively lower than that in the case of at- 
mospheric flow. Also, the total record length is usually 
much longer than the large-eddy turnover-time. In an 
atmospheric flow, on the other hand, the Reynolds num- 
ber is usually very large, but the characteristic time of 
the macroscale shear quite often becomes comparable to 
or longer than the record length. Thus, the homogene- 
ity and isotropy conditions are not well satisfied, and the 
problem of statistical convergence at large scales arises.™ 

Most studies of the scaling of structure functions in 
the DNS have been made using a very limited number 
of samples, and the results may therefore suffer from the 
effects discussed above. Hence, there are several factors 
to be checked before drawing definite conclusions from 
simulations about scaling; they are: (1) the degree of 
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isotropy of turbulence in the scaling region, (2) the ef- 
fect of external forcing, and (3) the effect of decreasing 
intensity of fluctuations. Forced simulation cannot avoid 
the first two factors, while the first and third factors 
are crucial in decaying simulation. Thus, we need to in- 
vestigate the effects of limiting the Reynolds number to 
finite values, and then determine the resulting longitudi- 
nal structure function. 

In the present paper we carefully examine the above 
three factors in decaying and forced simulated turbulence 
with various Reynolds numbers. In particular, flows of 
Reynolds number R\ = 70 are carefully studied. Twelve 
decaying turbulences with the same initial energy spec- 
trum (but different realizations) are simulated, and the 
forced turbulence is executed for a long period of time, 
providing an average over 126 samples, during approx- 
imately 50 eddy turnover times. The degree of local 
isotropy is studied through the incompressible-condition 
restriction on the second and third order structure func- 
tions. How external forces penetrate into the upper 
inertial region in forced turbulence and how the non- 
stationarity of the fluctuations deteriorates the upper 
inertial region in decaying turbulence are investigated 
through the Karman-Howarth-Kolmogoiov relation (re- 
ferred to as the KHK relation hereafter)© The KHK re- 
lation is an exact expression connecting the third-order 
structure function to the second order structure func- 
tion. On the basis of these studies, we estimate the 
scaling exponents of the longitudinal velocity structure 
functions in forced and decaying turbulence by^pxploy- 
ing the extended-self-similarity (ESS) method. li3EJ We 
find that the scaling exponents of the forced turbulence 
agree with other currently reported values.^ but that the 
scaling exponents of decaying turbulence are smaller for 
higher orders than the corresponding exponents of steady 
turbulence. A reason for this behavior is suggested using 
the equation for the structure function of arbitrary or- 
der, which is a generalization of the third-order structure 
function. 

Our paper is organized as follows. In section 2, the 
relevant data from several simulations of decaying and 
forced turbulence are summarized. Section 3 is devoted 
to discussion of the isotropy of turbulence, while the 
KHK relation is investigated in section 4. In section 
5 the longitudinal structure function is obtained in the 
scale range that is appropriate for the KHK relation and 
the conditions of isotropy to hold. The scaling exponents 
of longitudinal structure functions up to the eighth order 
are derived using the ESS method. Evidence that the ex- 
ponents for decaying turbulence are smaller than those 
for forced turbulence is also presented. In section 6, the 
effect of the finite Reynolds number on the structure of 
the local exponents is discussed, and a justification for 
the ESS method is given. In the appendices, the equation 
for the structure function of arbitrary order is derived 
and, in particular, the second-order equation is shown to 
reduce to the KHK relation. Additionally, the contribu- 
tion of the pressure gradient term in the cquatioa-for the 
structure function is shown to be short ranged.EJ> 





Run 1 


Run 2 


Run 3 


Run 4 


Run 5 


condition 


decay 


decay 


decay 


forced 


forced 


mesh points 


256 3 


256 3 


512 3 


256 3 


512 3 


Rx 


70 


90 


120 


70 


125 




1.11 


1.26 


1.42 


2.32 


2.03 


N 


12 


1 


1 


126 


45 



Table I. Various parameters in decaying and forced simulations. 
N is the number of samples over which the average is taken. 



§2. Numerical Procedure 

Three kinds of decaying runs were carried out for the 
Gaussian, random, initial-velocity field, using the same 
energy spectrum E(k, 0) = c(k/k ) 4 exp[— 2(k/k ) 2 ](ko = 
3 for Run 1 and ko = 1 for Run 2 and Run 3), but with 
different realizations. The number of mesh points used 
were 256 3 for Runs 1 and 2, and 512 3 for Run 3. Run 1 
consisted of 12 different sub-runs with the same Reynolds 
number. The numerical algorithm used for the calcula- 
tions is the pseudo-spectral method in space, and the 
fourth-order Runge-Kutta-Gill method in time. Aliasing 
errors were eliminated through application of the shift 
method. The computations were carried out on a par- 
allel vector- machine at NAL's Numerical Wind Tunnel 
and at RIKEN's Advanced Computing Center. The var- 
ious data were taken at the time where the energy dis- 
sipation had reached the maximum value. The relevant 
data are shown in Table I. The spatial resolution of the 
DNS's may be estimated by the value of the product 
fcmax?7, where rj is the Kolmogorov length. For all runs, 
the values obtained are larger than unity, implying that 
sufficient resolution was obtained for even the smallest 
scale. 

The numerical method used for forced turbulence is es- 
sentially the same as that used for decaying turbulence. 
The initial conditions applied to the forced-turbulence 
experiment were also the same as those used in the 
decaying-turbulence test. The random force is statisti- 
cally homogeneous, isotropic and Gaussian white, and its 
spectrum support is limited to the band 2<fc<3. Also, 
the spectrum form is constant in the band. The data 
are shown in Table I. After the transient period had 
passed, usually after a few turnover times, steady states 
were attained and the time average was computed over 
some multiple of the turnover time. For example, Run 4 
was averaged over 126 samples during 50 eddy turnover 
times, and Run 5 was averaged over 45 samples during 9 
eddy turnover times. The computations were performed 
at RIKEN's Advanced Computing Center. 

§3. Isotropy 

There are many ways to check the isotropy of turbu- 
lence. One is to examine the relation (u 2 ) — (t> 2 ) = 
(w 2 ), which reflects the isotropy of the largest scale, or 
the 'global" isotropy. In order to check the degree of local 
isotropy of eddies of scale r, which is necessary to com- 
pute the scaling exponents of the structure functions, 
we study the isotropic nature of the second and third 
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order structure functions. When the incompressible con- 
dition is applied to the second and third order structure 
functions, the following relations result for the isotropic 
case:EP 

rdD LL (r) 



D TT {r) = D LL (r) 



2 dr 



D LTT (r) = \ ^-rD LLL (r), 
6 dr 



(3.1) 
(3.2) 



where D LL = ((<5^) 2 ), D TT = {{Su t ) 2 ), D LLL = 
((<5wz) 3 ), Dltt = (Sui(5ut) 2 ) , Sui is the longitudinal 
velocity increment over a distance r, and dut is the 
transversal velocity increment. Instead of the differen- 
tial forms shown above, in the following derivations we 
employ the integral forms 

(D TT {r') - ^ LL (r')j dr' = ^rD LL (r), (3.3) 



r 1 

/ D LT T(r')dr' = -rD LLL (r), 
Jo 6 



(3.4) 



because numerical computation is less noisy for the inte- 
grals than for the derivatives. 

3.1 Decaying turbulence 

Figures 1(a) and 1(b) show comparisons with varying 
r for Run 1, averaged over 12 samples of different realiza- 
tions. This shows that the above relations are satisfied. 
The relation (3.4) for the third-order structure function 
is less accurate than that for the second-order, but it still 
holds up to r, = r/n = 160. (Note that from now on, 
the scale r will be expressed in units of the Kolmogorov 
length n.) The local isotropy of the structure functions 
is improved with an increase in the number of samples. 
To witness this, consider the isotropic relation for Run 
2, consisting of a single sample. In this case, the third- 
order relation was found to be slightly violated beyond 
r„ = 30 (figures not shown). 

Later we will determine the scaling indices of the lon- 
gitudinal structure functions in the inertial range us- 
ing the ESS method. We will employ the scale range 
16.0 < r* < 29.4 for Run 1. Because of this choice of 
scale range, the structure functions in the inertial region 
are not affected by any lack of accuracy in the degree 
of the isotropic condition of the turbulence. This is the 
case for Runs 2 and 3. 

Although it is not as easy to check the isotropy of the 
higher-order structure functions as it is to check the sec- 
ond and third-order functions, we expect the statistical 
convergence for the isotropy of the structure functions to 
be gradually lost as the order of the function increases. 
This is because the higher-order structure functions are 
dominated by rare, strong events that are oriented in 
specific directions. However, this does not necessarily 
mean that the degree of anisotropy of the fourth-order 
structure function is larger than that of the third-order 
structure function. The statistical convergence of even- 
ordered structure functions is faster than that of odd- 
ordered structure functions. 

In decaying turbulence, the state of isotropy is satisfac- 
torily established. Although a sufficient number of sam- 



ples is desirable for computation of the structure func- 
tions, the data of a single snap-shot (as in Run 2 or 3) 
are enough to determine the structure functions with re- 
spect to the degree of isotropy. This, however, is not the 
case with forced turbulence, as is discussed below. 
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Fig. 1. The isotropic conditions for Run 1, consisting of 12 iden- 
tical samples with different realizations, (a) The second-order 
structure function; (b) the third-order structure function. 



3.2 Forced turbulence 

It is important to note that the conditions (3.3) and 
(3.4) are not satisfied for forced turbulence as well as 
they are for decaying turbulence, particularly when the 
data of only one snap shot are processed. To see this 
explicitly, we show in Figs. 2(a) and 2(b) the curves 
of eqs. (3.3) and (3.4) for Run 4 at a certain instant 
after the turbulence reaches the steady state. Although 
the second-order structure function satisfies the isotropic 
condition, the same condition applied to the third-order 
structure function is violated at all scales. This indicates 
that an average over a prolonged period is needed to es- 
tablish the degree of isotropy for forced turbulence. On 
the other hand, Figs. 3(a) and 3(b) show curves for Run 
4 averaged over 126 samples, with the averaging-time be- 
ing approximately 50 eddy turn-over times. The isotropy 
condition for the second-order structure function is thus 
ensured quite satisfactorily. The isotropy condition for 
the third-order structure function, however, is only sat- 
isfied up to r* = 50. We expect that the condition for 
statistical isotropy becomes better satisfied as the num- 
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ber of samples increases, and as the averaging-time in- 
creases. This would suggest that steady turbulence with 
isotropic forcing needs to be studied over a sufficiently 
prolonged averaging time. 
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Fig. 3. The isotropic conditions for Run 4 averaged over 126 
samples during approximately 50 eddy-turnover times, (a) The 
second-order structure function; (b) the third-order structure 
function. 



Fig. 2. The isotropic conditions at one snap-shot for Run 4, when 
the energy dissipation reaches a constant level, (a) The second- 
order structure function; (b) the third-order structure function. 



§4. The Karman-Howarth-Kolmogorov Equa- 
tion 

In_this section we consider, through the KHK rela- 
tionJiiP how external forcing (in forced turbulence) and 
non-stationarity of the intensity of fluctuations (in de- 
caying turbulence) affect the transfer of energy in scale- 
space. The effects of the external forcing and non- 
stationarity on a weighted transfer of energy, which is 
the average of the product of the energy-dissipation rate 
and the multiples of the velocity-increment, associated 
with the higher-order structure functions are discussed 
later. 

The KHK relation is an exact equation, which is very 
rare in turbulence theory. The equation for the second- 
order structure function employing the isotropic condi- 
tion is written: 



D 



LLL 



(r) 



-sr - 



dD 



LL 



dr 



+ G(r)+F{r), (4.1) 



where 



and 



G{r) = -j 



9Dll ,] 

at 



dr', 



F(r) = —e m {k e r) 2 r, 



(4.2) 



(4.3) 



as is derived in eq. (A-9) of Appendix A. Note that 
equation (4.1) is valid for decaying turbulence as well as 
for the forced variety. Here e is the average dissipation 
rate, while Si n is the energy input rate due to external 
random forces, which are distributed about k ~ k e . Note 
that v represents the kinematical viscosity. The second 
term on the right-hand side (RHS) of eq. (4.1) expresses 
the effect of viscous damping, while the third term on the 
RHS relates the effect of the non-stationarity (decreasing 
in time) of the second-order structure function, which is 
significant for decaying turbulence. The fourth term on 
the RHS expresses the effect of the external forces. Note 
that only in steady, forced turbulence does e = £j„. It 
is of some interest to note that the first term on the 



RHS of eq. ([hi]) is negative, while the viscous term, the 
external-forcing term and the non-stationarity term are 
positive. 

In the incrtial region of stationary-turbulence with in- 
finite Reynolds number, one can neglect the second, third 
and fourth terms on the RHS of the above equation, so 
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that one gets Kolmogorov's well-known 4/5 lawH 

D LLL (r) = -\er. (4.4) 
5 

The inertial region in real flows with finite Reynolds 
numbers is regarded as having the relation (4.4) hold 
approximately, and hence the universal scaling property 
of the structure functions in the inertial region must be 
examined only under these conditions. In actual flows, 
however, the Reynolds number is not large enough to 
guarantee the wide spatial range in which eq. (4.4) will 
hold. 

For flows of finite Reynolds numbers, therefore, the 
important issue is how the transfer term — (4/5)er is af- 
fected by (1) the viscous damping term, (2) the non- 
stationarity term in decaying turbulence and (3) the ex- 
ternal forcing term in forced turbulence. Case (1) is sig- 
nificant in the lower inertial region, while (2) and (3) are 
pertinent in the upper inertial region. We investigate 
each term in the KHK equation for both the decaying- 
and forced-turbulence scenarios. The failure of eq. (4.1) 
is caused by the breakdown of the isotropic condition, 
and by the insufficient number of samples averaged. 

Now, concerning the numerical computation, e was 
calculated from the square of the spatial derivative of 
the velocity field, while the forcing term F(r) in eq. ( |4.3| ) 
was computed with the aid of the spectrum of the ran- 
dom forces, not from an external energy input at each 
instant. More precisely, we did not use the approxima- 
tion F(r) for the forcing term, but the exact expression 
(ATO) in Appendix A; we did this because, in this case, 
k e r is not limited to k e r -C 1. The other terms in the 
expression are all evaluated directly through processing 
the data. 

4-1 Decaying turbulence 

The Karman-Howarth-Kolmogorov relation now takes 
the form 

4 



D 



LLL 



(r) 



5 dr 
In Fig. 4, we depict — Dlll{t) /e and the co mbin ed sums 
of the terms on the right-hand side of eq. (4.5) divided 
by e, against the scale r* for Run 1 (averaged over 12 
samples.) The relation (4.5) is well satisfied for all values 
of r*. 

From Fig. 4 it can be seen that the peak of 
—DLLL{r)/(£r) is 0.53 at r* = 26.7, and the linear re- 
gion of Dlll(t) is extremely narrow around r* ~ 23. 
If the viscous effect is included, however, — Dlll(t-) + 
QvODll/ dr is in near-agreement with (4/5)er up to 

= 10. If the non-stationarity of Dll is taken into ac- 
count, agreement is guaranteed for any scale. The effect 
of the non-stationarity of Dll becomes more significant 
as the separation scale r increases, indicating that the 
amplitudes of larger-scale components decay faster than 
those of smaller-scale components. This is because the 
former do not have any energy supplied to them in the 
decaying turbulence. The non-stationarity term, which 
has been included in Fig. 4, is observed to decrease as 
r c for 5 < r* < 30, and where c is almost 3. This con- 



stant slope is surprising, because the scaling range of the 
non-stationarity term extends from the inertial to the 
dissipation region. 

We will now discuss the form of the KHK relation in 
connection with the structure function, with the details 
of the computations shown in the following section. As 
seen from Fig. 4, corresponding to Run 1, it is difficult 
to find the region where Dlll(t) = — (4/5)er. This 
implies that computation of the exponents in the inertial 
region by the direct-fitting method-is-not possible. We 
therefore employ the ESS methodO t3 to calculate the 
scaling exponents of the structure functions. The ESS 
method is known to extend .the scaling region down to 
the upper dissipation rangeJiU because the third-order 
structure function used in the abscissa contains the effect 
of the viscous-damping term. In the following section, on 
the other hand, we will show that the ESS method yields 
reasonable values beyond r* = 20 up to = 30, where, 
as Fig. 4 shows, the non-stationarity effect cannot be 
neglected. Thus the ESS method extends the scaling 
region upward in the scale-space in decaying turbulence. 

The KHK equation has been experimentally investi- 
gated in grid-turbulence by Danaila et a/.JIZP who took 
into account the non-stationarity term by approximating 
it using spatial differentiation under the Taylor hypoth- 
esis. They claimed that the KHK equation is correct up 
to ~ 200; however a slight discrepancy is seen around 
~ 20 (Fig. 3 in ref. 17). Such a deviation nega- 
tively influences the determination of the scaling expo- 
nents of the higher-order structure functions. The non- 
stationarity term takes the form r c with 1 < c < 2, in 
contrast to our assumed value of c = 3. The reason for 
this discrepancy is not currently known. 
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Fig. 4. The KHK relation for decaying turbulence for Run 1 av- 
eraged over 12 samples. 



4-2 Forced turbulence 

In this case the KHK equation takes the form 



D 



LLL 



(r) 



-er 



dD LL (r) 



F(r) 



(4.6) 



5 dr 
where the last term is numerically computed with the 
help of the exact expression (A- 10) in Appendix A. It is 
interesting to note that the relation (|4.6|) does not hold 
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when the data of only a single instantaneous snap-shot 
are processed. This situation is analogous to having the 
isotropic condition fail in Fig. 2. The reason for the sim- 
ilarity is that both are linked to a random feature in the 
external input. The last term in eq. flL6|) is calculated 
using an ensemble-average value from the spectrum of 
(|/i(fe)| 2 ), which is equivalent to a time average over a 
sufficiently long period. Let us now estimate how the in- 
put term varies from time to time. The number of input 
modes is N = 160 for Runs 4 and 5, so that a simple 
estimate of the relative amplitude of the input fluctua- 
tion is of the order 1/yN ~ 0.1. This implies that the 
external-forcing term at any given instant may have a 
relative error on the order of 10%. In order to obtain 
reliability of the forcing term to within a 1% error, we 
need 10 4 random forces, which corresponds to roughly 
one eddy-turnover time in our present simulation. To 



obtain the statistically meaningful result of eq. (4.6) we 
require a much longer averaging period than over only 
one eddy-turnover time. 

Figure 5 displays the KHK curve for Run 4 with its 
average over 126 samples, in a time interval of roughly 50 
eddy-turnover times. It is not surprising that the KHK 
relation holds here for almost any separation, since the 
relation is exact under the isotropic condition. Let us 
now inspect the contribution of each term in the KHK 
equation to its total value. According to Fig. 5, the peak 
of — Dlll(t) / (er) is 0.58 at = 26.8, and the flat re- 
gion of Dlll(t) jr is extremely narrow around r* ~ 26. 
This is very similar to the situation encountered in the 
decaying turbulence studied in Run 1. If the viscous ef- 
fect is also taken into account however, the agreement of 
—Dlll(t) + SvdDLL/dr with (4/5)er is nearly correct 
up to r* = 20. If the effect due to external forces is in- 
cluded, which begins to appear at r* ~ 20, the agreement 
between the two relations is perfect. 

Also of interest are the results of Run 5, averaged over 
45 samples, during nine eddy-turnover times. Here the 
KHK relation is not obeyed as well as in Run 4. Our com- 
putations confirm that agreement improves as the sample 
number increases, which explains why the agreement is 
better for Run 4. Note however that the improvement 
in agreement is small with modest increases in sample 
number. 

The KHK relation for a real flow was recently stud- 
icdli2P for low temperature helium gas flows, with R\ 
ranging from 120 to 1200, using an apparatus in which 
the turbulence is forced by two counter-rotating disks in 
a cylindrical container . Th e effects of large-scale motions 
were included in eq. (4.6) as the external-forcing term. 
The exact nature of eq. (4.6) was proved in this experi- 
ment, except for large scales, where non-isotropic forcing 
has a dominant effect. 



§5. Longitudinal Structure Functions 

Let us begin with an estimate of the maximum de- 
gree of the order of structure functions that we desire 
to measure accurately. We plot the probability-density 
function P{u r ) multiplied by it™ against u r , where u r is 
the longitudinal-velocity increment over the distance r. 
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Fig. 5. The KHK relation for Run 4 averaged over 126 samples. 
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As n increases, the peak of the bell shape of u™P(u r ) ap- 
proaches the larger value of |u r |. Beyond a certain value 
of n, the peak cannot be distinguished and rare strong 
events heavily influence the computation of (u"). Fig- 
ure 6 is an example of such a plot, showing uf.P(u r ) 
at i\ = 20.8 (Run 3). The abscissa in the plot is 
u r j \J 5*2 (r), while the ordinate is pre-divided by Si(r)^ . 
The curves decay quickly for large amplitudes. For the 
plots of u"P(u r ) with n > 9, we observed that the decay 
of the curves is not sufficiently fast to ensure that mo- 
ments higher than the ninth order are computed accu- 
rately (figure not shown). Thus we conclude from these 
simulations that n = 8 is the maximum order allowed for 
the condition of a sufficient-number-of-events to be met 
in the statistical sense. 

In this paper we will concentrate our efforts only on 
the study of the longitudinal- velocity structure function 



S n (r) 



(5.1) 



Before going on to determine the exponents of S n (r), we 
will examine the structure functions of odd orders. For 
this purpose, let us recall the probability-density func- 
tion (PDF) of u r : P(u r ). The peak of the PDF is located 
at the small positive value of u r , while its tail is more 
populated on the negative side than on the positive side, 
in accordance with (u r ) = J u r P(u r )du r — 0. For (u^), 
on the other hand, the contribution from the negative 
values of u r is larger than that from the positive values, 
implying a negative value for S3 (r). Note that S n (r) 
becomes more negative as the odd values of n increase. 
Since the structure function of odd order is determined 
from the difference of two terms of almost equal mag- 
nitude, contributed from positive and negative values of 
delicate balance between the two terms is responsi- 
ble for the behavior of the odd-order structure functions. 
The structure functions of even order are the sum of 
those same terms, and do not rely on that intricate bal- 
ance. Thus reliable evaluation of the odd-ordered struc- 
ture functions is more difficult than evaluation of the 
even-ordered structure functions, and a smooth change 
in the exponents of S n (r) with respect to varying r in 
n, is expected only for flows of extremely large Reynolds 
number. A way of overcoming this difficulty for moderate 
Reynolds numbers is to employ the generalized structure 
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functions (|Mr| n 
thus introduce 



instead of S n (r) 



S„(r) = (\u r 



We 



(5.2) 



S n (r) will now reflect the change of the entire form of 
PDF, with respect to u r , so that the scaling exponents 
of S„(r) are now expected to vary smoothly with n. It 
is expected that the exponents of S n (r) and S n (r) for 
an odd number of n will be identical for large Reynolds 
numbers. In what follows, we drop the symbol tilde on 
S n so that no confusion will result. 

When we plot S n (r) against r*, it is difficult to read 
the exponent £„ defined in the relationship 



S n (r) 



(5.3) 



An example of this is shown in Fig. 7, where curves 
of 5*6 (r) for Run 1 (decaying turbulence, 12 samples) 
and for Run 4 (forced turbulence, 126 samples) are plot- 
ted against r* . The magnitude of both these quantities 
changes in such a way that both curves agree with each 
other in the dissipative region. It should be noted that 
the slope of S 6 (r) for decaying turbulence (Run 1) is 
smaller than that of the same quantity for forced turbu- 
lence (Run 4). The local slope dlogS 6 (r)/dlogr never 
assumes a constant value over an extended range, as 
shown in Fig. 8. Often- for moderate Reynolds number 
flows, the ESS methocO'EJ* is used instead of the above 
llpwing the work of many authors 
J0> we have plotted S n (r) = (\u r \ n ) 
against S^r) = (|M r | 3 ) instead of (m 3 ), for the reason 
stated previously in this section. Figure 9 shows a plot 
of Se(r) against the S3 (r) corresponding to Fig. 7; the 
power law tendency is clearly seen in this comparison. 

The exponents of S n (r) are obtained from the graph of 
the local exponent of log S„(r) against log S3 (r). Com- 
bining S„ (r) ~ r^ n and S3 (r) ~ r^ 3 , we have 




S„(r) ~ S 8 (r)<»/& 



(5.4) 



Hence, the slope of the plot of log S„(r) against log Ss(r) 
is Cn/C3 = Cnj since C3 is assumed to be unity. 
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Fig. 7. The plot of Se(r) against r t for Run 1 (decaying turbu- 
lence) and Run 4 (forced turbulence). 
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Fig. 8. The plot of the local slopes of S${r) and Se (r) against r* 
for Run 1 (decaying turbulence) and Run 4 (forced turbulence). 
The power-law dependence is not appreciable. 




Fig. 6. The plot of U%P( u r ) against U r at r, = 20.8 for Run 2, 
where U r = u r / ^Js^ir). 
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Fig. 9. The plot of Se(r) against Si(r) for Run 1 (decaying tur- 
bulence) and Run 4 (forced turbulence). Here the power-law 
dependence is appreciable. 
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5.1 Decaying turbulence 

In this section, we will employ the local slope method, 
which has the advantage that it requires less care in 
choosing where to draw the appropriate straight fitting- 
line on the graph. In Fig. 10 the local slope of log S n (r) 
of Run 1 is depicted. Note that the abscissa in this plot 
is not log S3 (r), but logr*, which is determined from 
Ss(r). It is remarkable that all the local slopes take on 
constant values in the interval 16.0 < r* < 29.4. This in- 
dicates that the scaling determined here is reliable. Runs 
2 and 3, using different Reynolds numbers and different 
realizations, have almost the same exponents as Run 1, 
although the scaling-range differs slightly from case to 
case. The local slopes £„ for Runs 1 to 3 are listed in 
Table II along with the scaling ranges. The last entry 
in the tablfi shows the standard numbers fitted by an 
SL model.Er Observe that when the Reynolds number 
increases, (2 increases while ( n , (n > 4) decreases. The 
measured exponents in the decaying turbulence simula- 
tions are slightly (but definitely) less than the standard 
numbers obtained in forced turbulence studies. 

We will now turn to a discussion of the condition of 
isotropy and the extent of non-stationarity in connec- 
tion with the determination of the exponents of struc- 
ture function. The isotropic conditions for the second- 
and third-order structure functions are guaranteed in the 
scaling range. We will now seek a connection with the 
non-stationarity of the functions. Such effects are in- 
ferred from the study of the KHK equation, or equiv- 
alently, from the third-order structure function. The 
KHK relation of Run 1 (Fig. 4) shows that the effect 
of the non-stationarity term dDr,L(r)/dt is not negligi- 
ble in the scaling range 16.0 < r* < 29.4, which is the 
range in which the scaling exponents of the structure 
functions for Run 1 are satisfactorily obtained. Study of 
the results reveals that the non-stationary effect is 22.4% 
of the (4/5)er at r* = 29.4, which is the upper cutoff, 
while it is only 8.2% at r* = 16.0. If the non-stationary 
effect is included in the total interaction, the energy flux 
in the scale-space of r is 



e(r) = e 1 



15 f r dD LL {r',t) 



Aer 5 



dt 



4 dr' 



(5.5) 



which can be derived from eq. (4.5). Since the integrand 
of the second term on the right-hand side of the above 
equation is negative in decaying turbulence, it follows 
that the function s(r) decreases as r increases. 

In order to see the effect of non-stationarity on higher- 
order structure functions, it is useful to examine the 
equation for the n + 1th moment of the increment of 
the ith velocity component Wi = Ui{x2) — Ui(xi), which 
is derived in Appendix A as eq. (A- 7). It is: 



d_ 

dr k 



- — (wkW™) + n( w. 



G^{r)-D^{r) + 2vVl( 



where G$ represents the non-stationarity effect 



(5.6) 



(5.7) 



and D$ expresses the correlation of the dissipation rate 
and u>", which is written: 



DW(r) = 2n(n-l)( £i (X,r)wr 2 ). 



(5.8) 



Here is the dissipation rate of the kinetic energy of the 
ith velocity component at the two points X ± r/2. It is 
written: 

e i (X,r) = J[|V x u; i | 2 + 4|V r ^| 2 ] 

= V - [\Vui(X + r/2)| 2 + \Vu % (X - r/2)| 2 ] . (5.9) 
A weighted energy flux can then be written as 

< £l (X,rK"- 2 )[l-pW(r)], (5.10) 

where 

1 



2n{n~l)( £l {X,r)wr 2 y 



(5.11) 



Since Ei(X,r) is positive-definite, as can be established 
from its definition (5.9), and the PDF for Wi is negatively 
skewed, then the denominator in p„ (r) is positive for 
even n and negative for odd n. On the other hand, the 
numerator is negative for even n, and positive for odd 
n, because its amplitude decreases with time for decay- 
ing turbulence. Thus p„ (r) is a positive quantity. The 
weighted flux is thus smaller than the value without the 
inclusion of the non-stationarity effect. Therefore it is 
possible that the non-stationarity effect yields a differ- 
ent scaling for the structure functions in question. 

Although the KHK equation was not computed for 
Runs 2 and 3 (because the data for the time deriva- 
tive were not stored during measurement), we can es- 
timate the effect of non-stationarity by comparing the 
term Dlll — 'ovdD^jdr to — (4/5)er, since the effect 
of decay is attributed to their difference. To save on 
space in this paper, we do not include this figure in the 
text, but we have noticed that the non-stationarity effect 
is not negligible in Runs 2 and 3 where larger Reynolds- 
number values arc exhibited. 




Fig. 10. The local slopes of the structure functions for Run 1; 
pluses stand for £2, crosses for £4, stars for £5, squares for ^6, 
triangles for £7, and circles for fg. The straight dotted lines are 
drawn with numbers averaged over the scaling region. 
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Run 1 


Run 2 


Run 3 


Run 4 


Run 5 


SL 


condition 


decay 


decay 


decay 


forced 


forced 




fitting region r* 


16.0 ~ 29.4 


11.7 ~ 26.0 


12.4 ~ 48.0 


11.5 ~ 28.1 


10.2 ~ 30.6 




(2 


0.703 


0.705 


0.713 


0.690 


0.692 


0.696 


C4 


1.266 


1.264 


1.252 


1.288 


1.284 


1.280 


Cs 


1.507 


1.502 


1.477 


1.555 


1.546 


1.548 


Ce 


1.724 


1.718 


1.680 


1.804 


1.788 


1.778 


C7 


1.920 


1.915 


1.864 


2.037 


2.011 


2.001 


Cs 


2.096 


2.095 


2.035 


2.254 


2.217 


2.211 



Table II. Scaling Exponents. The last entry, denoted as SL, shows the standard values presented in ref. 3. 



5.2 Forced turbulence 

We will now discuss the structure functions for Run 
4 averaged over 126 samples for approximately 50 eddy- 
turnover times. The calculated local slopes have been 
plotted against r* and are shown in Fig. 11. The con- 
stancy of the slopes in the interval 11.5 < r* < 28.1 is 
extremely impressive. Although the Reynolds number is 
70, the exponents are very close to the standard values. 
Run 5, which has a flow of even higher Reynolds number, 
yields similar results, as can be seen in Table II. Note 
that the scaling range is rather narrow in Run 5, despite 
its larger Reynolds number, than in Run 4. We suspect 
that the reason for this is that the sample number is not 
large enough in Run 5, where we averaged over only 9 
eddy turnover times, compared to the 50 eddy turnover 
times used in Run 4. In accord with current understand- 
ing of the process, we see that longer runs yield wider 
scaling ranges. 

We will now investigate the connection between 
isotropy and external forcing. Note that for the scal- 
ing region used for this test, the isotropic conditions 
are well satisfied, while the KHK equation is affected 
slightly by the external- forcing terms. According to Fig. 
5, the effect of external forces is small in the above scal- 
ing range. However, careful inspection of the numerical 
results shows that the forcing term in KHK is 13.3% of 
(4/5)er at r* = 28.1 and 2.3% at r* = 11.5. This reveals 
that external forcing does not affect the lower inertial 
region, which is where the constant local slopes begin, 
but some effect is evident in the upper inertial region. 

In forced turbulence, the equation for the n + 1th mo- 
ment of the increment takes the form: 

= H^\r)-D < i\r) + 2vVUO, (5-12) 

where ifj, represents the contribution of the external 
force, and is written as 

ffW(r) = l -n{n - l) £m (fc e r) 2 R~ 2 )- (5.13) 

(See Appendix A for a more complete explanation of this 
term.) The importance of the external forces can then 
be estimated through the ratio of Hn \r) to the transfer 
term on the left-hand side of eq. (5.12): 



where a n = 1 + Cn-2 ~ Cn+i, is an increasing function of 
n and approaches the constant value 2/3, according to 
the SL model. This estimate predicts that the influence 
of the external forces will diminish as the order of the 
structure function increases. Consequently, the effect of 
the external forces will not be significant in the higher- 
order structure functions. 




10° 10 1 10 ; 

n 



Fig. 11. The local slopes of the structure functions for Run 4. 
The plus symbols represent £2, the crosses £4, the stars £5, the 
squares fg, the triangles £7, and the circles £g. The straight 
dotted lines are drawn with numbers averaged over the scaling 
region. 



5. 3 Comparison 

The results of our careful simulations show the scaling 
exponents to be definitely smaller in the decaying tur- 
bulence simulations than in the forced turbulence simu- 
lations. Recently, Boratav and PelaP obtained the scal- 
ing indices for simulation of decaying turbulence with 
R\ = 82, which agrees with the standard values. How- 
ever, they did not employ a method of computing the 
local slope of the structure functions and identifying the 
scaling region with the constant local slope. They simply 
drew a best-fit line on their graphed curve to estimate 
its slope. Such a method will predict a different set of 
scaling indices depending on which region is chosen. 

Let us now consider the reasons why the exponents 
are smaller for decaying turbulence than for forced tur- 
bulence. One conceivable reason for this phenomenon 
is that it is caused by the effect of the system's non- 
stationarity, as discussed in §5.1. One may argue that 
large Reynolds number flows do not suffer from non- 
stationarity effects; however, since the magnitude of the 
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non-stationarity effects decreases only by a certain power 
law as r decreases, the entire inertial region should be af- 
fected by the influence of non-stationarity. 

Another possible reason for this discrepancy is that 
the dissipative structure may be different in decaying 
turbulence than in forced turbulence. Quantitatively, 



< £l (x,r>r 2 ) 



(5.14) 



may scale differently in decaying- and forced-turbulence 
environments. As is shown in Appendix B, the pressure 
gradient term in eq. (5.6) is not affected by fluctuations 
at a far distance, and it obeys the same scaling relation- 
ship as the inertial term. Hence the left-hand side of eq. 
(5.6) is represented primarily by the inertial term. The 
essential equation for the inertial range is 

(1 + = -2n(n-l)( £l (J,rK- 2 ), 

w ^ . (5 - 15) 

where Il-i represents the contribution from the pressure 
gradient (see Appendix B). If the fluctuation of the dis- 
sipation rate is neglected, then the above equation can 
be replaced by: 



fl + jW ) — ( 
[L + in - l) dr k { 



-~n(n — l)s(w r i 



2 ), (5.16) 



which governs the K41 scaling. In real turbulence how- 
ever, the dissipation rate fluctuates in such a way that 
[uii(X , r)] n ~ 2 a.nd £i(X ,r) are positively correlated. Un- 
der these circumstances we have intermittent scaling, 
where the scaling exponents are smaller than their K41 
values. 

It is possible that the correlation of [wi(X, r)] n ~ 2 and 
Ei(X, r) differs for decaying turbulence and forced turbu- 
lence. This may happen because the dissipative structure 
differs for the two turbulences and t he n on-stationarity 
effect may also affect the correlation ( 5.14 ) . Precise anal- 
ysis of the smaller exponents in decaying turbulence is a 
topic left for future study. 

§6. Finite Reynolds Number Effect on the Scal- 
ing Exponents 

So far, we have computed the structure functions for 
decaying and forced turbulence using a statistically suffi- 
cient number of samples, and have estimated the scaling 
exponents using the ESS method. In this section, we dis- 
cuss the behavior of the structure function. If one looks 
at the plot of the structure function against log r* (as for 
log Se (r) in Fig. 7, for example,) it can be seen that the 
curve's shape may be approximated in the scaling region 
by a quadratic form of log : 



C 

log S n (r) = A n + B„logr* - -y (logr*) 2 . 



.1 



Here A n , B n and C n are positive. Note that as the 
Reynolds number, R\ 1 increases, the linear range of the 
above equation increases, indicating that C n decreases 
with R\. From the results of our studies, we suspect 
that B n varies little with changes in R\. The local ex- 
ponents can then be computed from 



As r* increases, the structure functions saturate (that 
is, become independent of r), as can be seen from the plot 
of Sq(t) in Fig. 7. This implies that Cn( r ) g° es to zero 
at the saturation scale. We then require that £„(?") = 
at r = Lq. Note that we do not necessarily think that 
I/o is the same as the external length. Consequently, 



L B n 
log — = — — = const 



M, 



and substituting eq. (|6.3D into eq. (S.2) yields 



Cn{r) =B n (l 



log 

M 



(6.3) 



(6.4) 



Since C n decreases with R\, M will increase with R\. 
Emp loying the ESS expression £„ = (r)/^3 (r) , eq. 
(6.4) leads us to 

g n (l-logr./A/) = Bn 
£ 3 (l-logr,/M) ~~ B 3 



= ^, (6-5) 



which is independent of r*, so long as we stay in the 
scaling range. This explains why the ESS works so well 
for our simulations. In the limit of infinite R\, B% — > 1 
and therefore, B n — > Q n . 

We will now discuss the evidence supporting eq. (3.4), 



relying on our present simulations and other experiments 
employing large Reynolds numbers. Figure 12 shows a 
plot of the local slopes Cn( r ) against for Run 4, in 
which various straight lines fitted in the inertial region 
(where the scaling exponents are to be determined) are 
extrapolated outside of the scaling region. All the lines 
converge to a single point, (10 M ,£„ = 0), although the 
data points for large r* values are su ppre ssed to empha- 
size the convergence. The relation (|6.4| ) is remarkably 
well satisfied with M = log 143 = 2.16. Other runs 
with forced turbulence show the same tendency with 
M = log 188 = 2.27 for Run 5. It is clear from these 
results that M increases as R\ increases. For the decay- 
ing turbulence of Run 1 using the same R\ as Run 4, 
M = log 112 = 2.05 which is slightly smaller than the 
corresponding value for Run 4. This is interpreted as 
likely being due to a possibly larger C n value for decay- 
ing turbulence resulting from the non-stationarity effect. 

The linearity of £„ (r) with respect to log r* is demon- 
strated in experiments for turbulent flow£_wi±h Reynolds 
numbers larger than the present studyllS'EP It should 
be noted though that the scale ranges where the linear- 
ity holds in these experiments are located at larger scale 
values than in our studies. Note that the local slopes 
are shown for more restricted values of n in those other 
experiments; for example: M = log 2 x 10 5 = 5.28 in the 
flow of R\ = 2000,E2P where M is read from the plot of 
Ce(r) (see Fig.^(d) in ref. 20). According to the flow 
of i? A = 195000) M ~ 11 is estimated from ( 2 (r) (horn 
Fig. 4(a) in ref. 16), and M ~ 10 from ^(r) (from Fig. 
4(b) in the same publication). In order to see how M 
varies with R\, we create a log-log plot of M against 
R\, as in Fig. 13. From this we can establish the rough 
relationship: M ~ R\' 3 - We expect this relationship to 
be confirmed in future work. 



C»(r) =B n - C„logr* 



(6.2) 
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Fig. 12. The verification of relation ( ]6.4j ) for Run 4. The data 
points for large values of r, are suppressed to emphasize that 
the scaling lines converge to a single point at (r» = 143, C, n = 0). 
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Fig. 13. The curve of M against Rx- A circle symbol represents 
the data for Run 4, and a square the data for Run &ri A plus 
indicates the data for the experiment of R\ = 2000£lP and a 
cross (derived from £2(0) and a triangle (derived from C&( r )) 
show the results for experiments of R\ = 19500.|i3^ The broken 
line represents the relation M ~ 



§7. Summary and Discussion 

In this paper we have examined turbulence condi- 
tions to accurately determine the longitudinal structure 
function through the incompressible expressions of the 
second- and third-order structure functions and through 
the KHK relation. The examination of the isotropy 
for the second- and third-order structure functions, in 
terms of the incompressibility conditions, showed that 
the isotropy of these functions is well satisfied when av- 
erages are taken over the ensemble of the initial field 
or over a sufficient period. It was found that the non- 
stationarity effect due to the decaying and forcing effects 
penetrates into the upper inertial range, and the ESS 
method must be applied to measure the scaling expo- 
nents of the structure functions on scales smaller than 
those of penetration. 

The degree of anisotropy of turbulent flow will increase 
with the order of the structure function, as can be in- 
ferred from our observation that the third-order struc- 
ture function is more anisotropic than the second-order 
function. The problem lies in determining to what extent 



the anisotropy affects the higher-order structure func- 
tions. It is difficult to numerically predict the anisotropic 
effects, because the simple expressions for the isotropy of 
such systems, as related in eqs. (3.1) and (3.2), have not 
yet been extended to include the higher-order structure 
functions. Also, we do not yet know how the effect of 
non-stationarity changes with the order of the structure 
functions. This could possibly be estimated by examin- 
ing numerically the factor defined in eq. ( 5.11 ). 

We have computed the scaling indices of the structure 
functions up to the eighth order, which we believe to be 
the maximum order for which accurate results are given 
under the above considerations. It was found that the 
exponents are definitely smaller for decaying turbulence 
than for forced turbulence. Further studies are necessary, 
however, to derive definite conclusions either extending 
or disproving this observation for turbulence with higher- 
order structure functions. Our limited current results for 
simulated decaying turbulence, however, do indeed show 
the smaller numbers. 

Since smaller exponents correspond to the more in- 
termittent characteristics in the turbulence, we see that 
decaying turbulence is more intermittent than forced tur- 
bulence. It is widely recognized that a coherent struc- 
ture with a very large velocity gradient is of the form 
of a vortical tube which has a diameter and length of 
the order of the Kolmogorov- and integral- lengths. It 
is plausible that random forcing at a large scale acts to 
destroy the coherence of such a structure over the inte- 
gral scale. On the other hand, in decaying turbulence 
there is no such mechanism to prevent the development 
of a coherent structure in the velocity field, other than 
basic restrictions imposed by the computational box-size 
or experimental geometry. For this reason, we observe 
that decaying turbulence is more intermittent than the 
forced variety. 

We have argued that the local-scaling exponents ( n (r) 
can be expressed as a linear function of logr*, reflect- 
ing the finite nature of the Reynolds number. When 
the Reynolds number becomes large, the correction term 
arising from the effect of the finite Reynolds number van- 
ishes. 
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Appendix A: Derivation of Karman-Howarth- 
Kolmogorov Equation 

Although the derivation of the KHK equation is-given 
in most standard books dealing with turbulenceJIIP we 
will derive this equation for a structure function of arbi- 
trary order and will include the effects of forcing. Note 
that the KHK equation turns out to be nothing more 
than the equation for the second-order structure func- 
tion. 

Let Ui and «2 be the velocity of the flow at locations 
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Xi and x 2 . The average velocity V and the relative 
velocity w are defined as 

w(X2,Xi) = M 2 — Ml, V = (u 2 + Ui)/2. 

The equation for w is derived by taking the difference of 
the Navier-Stokes equations at x 2 and x\\ 

d 

-w(x 2 ,x 1 ) + (u 2 ■ V 2 + mi • Vi)i»(jc2, oji) 



dt 



d : _d_ 

dx 2 dx\ 



(p(x 2 ) -p(xi)) 



+z/(V 2 + v;)to(aj 2 ,xi) + f{x 2 ) - f( Xl ), (A-l) 

where / is an external random Gaussian force, p is the 
pressure and v is the kinematical viscosity. Notice that 
the variable x 2 is held constant when the partial deriva- 
tive with respect to x\ is taken, and similarly for x\ 
when the partial derivative with respect to x 2 is taken. 
Now, we introduce the quantities 



x 2 - X\. 



1 d d_ 
2dX + dr' 



X = (aji + aj 2 )/2, r: 

Since 

did d d 
dxi 2 dX dr' dx 2 
we are led to the expressions: 

m 2 • V 2 + mi • Vi = V(X, r, t)-V x + w(X, r, t) • V r , 
V 2 + V! = v x , 
1. 



v 2 , 



-vV 

2 v 



2V£. 



Here r) is the velocity difference w(x\, x 2 ) with a 

change of variables from x\, x 2 to X, r. 

On substituting the above relations into eq. (A-l), 
the Navier-Stokes equation for the velocity difference be- 
comes 



Wi(X,r,t) 



d d 
= -w k (X, r, t)—Wi(X, r, t) - — <5p(X, r, t) 



2 A 



2V 2 U i (X,r,i)+5/ i (X,r,i), 



(A2) 



where 

<5p(X, r, i) = p(jc 2 ,*) - p(xi,t), 
5f(X,r,t) = f(x 2 ,t)-f(x u t). 
Note that the following incompressible condition 

ik- y = \ik- {u{X2)+u{xl)) 



i 



d 



d 



9a; i 



0. 



holds, because r is held constant when the partial deriva- 
tive with respect to X is taken. Similarly, 



d d Tr d 

• w — — • V = — • w 

dX dr dr 



0. 



Multiplying eq. (A-2) by w™ 1 , and reorganizing the 
dissipation term slightly, we get 

I| < + Iv(X,r,i).V x < 



— w ■ V r w" 
n 



"2 v x • « 



n-l 



^^(X,r,i) 



9X 

n — 2 , r <* „..n— 1 



-iA7 2 t(;, n 
n 



Here, 



-2(n- l)Ei{X,r)w?- 2 + 8f iW . 



Si{X,r) = - [\V XWl \ 2 +4\V r wi\ 



(A3) 



= - [\V Ui (X + r/2)| 2 + \V Ui (X - r/2)| 2 ] (A4) 

is the summation of the energy-dissipation rate of the 
ith component of the velocity field at the two points x\ 
and x 2 . 

Let us now calculate the correlation of Sfc and w"" 1 , 
which appears in eq. (A-3), 

F i (n) = (6f i (t)wr 1 ) f , 

where {•) j signifies the average over Gaussian random 
forces. Since the random force is simply delta correlated 
in time, and Sfi(s)(s < t) correlates with 5fi(t) only at 
t = s, it is appropriate to use the approximation 

Wi(t) = / 6fi(s)ds 

J — oo 

in the computation of the correlation (5fi{i)w™~ ls ) ^ 
Substituting this result into one of the n—1 components 
of Wi, we have 

Fi(n) = (n - l)wr 2 f (SMt)5M S )) f ds. 

J — oo 

In the present simulation, external forces are added 
randomly, once per time-step At, so that 

F l (n) = (n-l)wr 2 (8.fi) f . 

The correlation of random forces is calculated 
through 

<*/?>, - E (Mm(Q)) f ^ k+q) - x 

k.q 

x (^k r/ 2 _ c -tk-r/2^ _ & -iq-r/2^ 



= E^( fe )/*(fe))/(l-e- <fe - r )(l 
k 

^2^(|/ 4 (fe)| 2 ) / (l-cosfc.r), 
k 

where the homogeneity of the system is employed. Since 
the random forces are limited to wavenumbcrs smaller 
than 1 /r, where r is our scale of interest, cos k ■ r may be 
expanded in terms of kr. Furthermore, we can make use 
of the property that the external forces are distributed 
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in an isotropic way. Then 



(Sff) f = E (\m\ 2 ) f {k ■rf='-Y J (\Mk)\ 2 ) f k 2 

k 3 k 

When wc introduce the definition 

£ fc <l/*(fc)IV 2 



k 2 = 



we have 



{5ff) f = \{k e rYY.(\fm 2 ) r 



k 

Since the energy-input rate £,•„ due to the external forces 
is 



e m = 3(f?) f = 3j2(\Mk)\ 2 ) 



we are led to the relation 
1 



F i( n ) = n(n - l)e m (k e rY 



(A-5) 



(A-6) 



Substituting eq. ( A-6 ) into eq. (A-3) and taking the 
ensemble average over the velocity field, we have 



d i n\ d i 

(W, ) = [WkW 



d 



dV 11 dr k 

+2 V V 2 r {w™)-2n{n-l)(e l {X,r)w™- 2 ) 

+ l -n{n~l)e m {ker?(wr 2 )- (A- 7) 

Note that we have not made the summation over i at 
this point. 

Let us now focus on the KHK relation, which is derived 
by substituting n = 2 into eq. (A- 7) and making the 
summation over i. We then have: 

3 E(^) + ^E(w 2 > 



dt 



= 2vV 2 r ]T (w 2 ) - 4]T {e t (X, r)) + -e m (k e rf 



(A- 



Making use of the relations from 







i 

i 



we arrive at 
J__9 3 d_ 

r 2 dr dt 



r 3 ^D LL 



J_d_ 

3r 2 dr 



ld_ 

r dr 



r A D 



LLL 



2v d 
r 2 dr 



i d 4 d 

r dr dr 



-4e+ -e m (k e rf 



First, we multiply the above equation by 3r 2 , and then 
integrate it over r from to r. Second, we multiply the 



result by r, and integrate it again over r. We then get 



Dlll = -T er 
5 



dD LL 



dDLL -r'\\r< + l £m (ker) 2 r, 



dt 



35 



(A-9) 



for (k e r) 2 <C 1. Note that the last term on the right-hand 
side of (A-9), for arbitrary r, is given by: 

E k tfm 2 )^ mk)n 
3 , 9 

x| — h —r — — sin kr + — ; — 7 cos kr 



k 3 r 3 



fc 4 7 



fc 5 r 5 



sin kr 



(A-10) 



Appendix B: Contribution of the Pressure Gra- 
dient Term 

The equation for the structure function of arbitrary- 
order has been derived above. Each term appearing in 
the equation acts locally in space except for the pressure- 
gradient term. The pressure is usually considered to be 
a long-ranged influence, because it is related to the far 
velocity-fields through the Poisson kernel. Consequently, 
the same is commonly thought to be true for pressure- 
gradient forces. Here we show that the pressure-gradient 
term is not as long-ranged as thought. The fact that the 
pressure-gradient term has a local nature was mentioned 
previously by L'vov and Procacciajij' but the argument 
can be better illustrated if we use the results in the Ap- 
pendix A. 

Let us take the divergence of eq. (A-2) with respect to 
Xi. Since the external forces are divergence-free, we are 
led to the equation: 



W 2 x Sp(X,r,t) 



d 
dX~ 



d 

w k (X,r,t)— Wj(X,r,t) 
dr k 







+V k (X,r,t)— Wj (X,r,t) 



(B-l) 



The source term can be simplified as follows. The in- 
comprcssibility condition reduces the source term to 

dw k (X,r,t) dw 3 {X,r,t) dVj(X,r,t) dw k (X,r,t) 



8X 4 



dr k 



dx k 



8X< 



where the interchanging of (J «-> k) has been done in the 
second term. Since 



d 
dXk~ 



Vj(X,r,t) 



1 d 



2dX k 
d 
dr k 

d_ 

dr k 

the source term becomes 

dw k (X,r,t) 8wj(X,r,t) 



MX + rV2) + U ,(X-r/2)] 



[ Uj (X + r/2)- Uj (X-r/2)] 



—Wj(X,r,t), 



dXj 



= -2- 



dr k 
_d d_ 

' dXj dr k 



w k (X,r,t)wj{X,r,t). 
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The pressure-gradient term then reduces to: 
JL Sp (X,r) = J dX'K^X.X 1 ) 



haves as (R/r)~ s , for R/r > 1 with S > 0. Thus the 
contribution from the upper boundary of the second in- 
tegral of eq. (B-7) can be neglected, because 



x— Wfc (X>,iH(JT,r,i), (B-2) 



dR 



2ttR 3 



, R i R 3 r 



RT 







as i? 



where 



K lj (X,X') = 



i a a 



i 



2tt dXidXj IX - x'\ 



(B-3) 



The pressure-gradient term appearing in eq. (A- 7) be- 
comes 

JJ„_ 1 ( i »,(r))E I! /,r 1 ^l' 



(B-9) 

Since this integral is dominated by scales of R compa- 
rable to r, s£ n_1 ^(i2 = r,r) can be approximated as a 
homogeneous function of r, 



,J dRK ia (R)B^- l) {R,r,t), (B-4) 



B<£-V{R,r,t) 



(X,r,t)-£-Wf,{X + R,r,t) 



drp 

xw a (X + R,r,t)\, 



(B-5) 



where the summation is taken only over the repeated 
Greek variables, and Kij(R) is a dipole-type interaction: 



K l3 {R) 



1 



2irR 3 



(B-6) 



The average value (B-5) does not contain the following 
contribution, 

(w^-\X,r,t))l^-wp{X + R,r,t)w a (X + R,r,t) 

since the latter average vanishes due to the fact that the 
divergence of the second-order structure function is zero. 
Therefore, the two quantities at X and X + R, that is: 

d 



^(ri/.r)*^- 1 ^!/,/*), 

where v = R/R, /x = r/r and (y, fl) is a geomet- 

ric function of v and fi. <j> is the same as the scaling 
exponent of (w^' 1 ^wpw a ^. 

Now E n -i(wi(r)) can be evaluated as 

£ n -iMr)) « r*nji n-1) (M), (B-10) 
J *"~ 1)(/X) = 2^/ dr!(iy) (3 ^ Q " <y »)4 n-1) 

(B-ll) 

This development implies that the greatest contribution 
to the integrals (B-4) comes from the range of R com- 
parable to r, meaning that the dominant contribution 
to the pressure-gradient term is from the near correla- 
tion on the order of R ~ r. As far as the scaling of the 
pressure-gradient term in r is concerned, the contribu- 
tion from the pressure-gradient in eq. (A-7) obeys the 
same scaling as the inertial term -^-{w a wf). We may 
then write: 



£„_iMr))=4^(M)^K<>, 



(B-12) 



where I„_i is a function of the order one, depending on 
fi. The factor n on the right hand side of eq. (B-10) is 



(X,r,t) and — — wa(X+R, r, t)w a (X+R, r, t) absorbed into the exponent of jf—(w a w, 
drp ar " 

must be correlated. It is certain that such a correlation 

function must decrease with increasing separation R for 
R > r. 

For \r\ in the inertial range, the integral of eq. (B-4) 
may be written as 

f dR+ f dR) K ia (R)B^ l - 1 \R,r). (B-7) 

JR<r JR>r / 

The first integral, B^ l ~ 1 \R 1 r), can be Taylor-expanded 
in \R/r\ as 

(R, r) = B^- 1 ) (0, r) + V ' R B^ (0, r) ■ R+ ■ ■ ■ . 

(B-8) 

Whe n th is form is substituted into the first integral of 
eq. (B-7), the first term of eq. (B_8) vanishes due to the 
symmetry of the system: 



dR 



1 



2ttR 3 



6 R? 6r - 



= 0. 



The other terms will make their contributions mostly 
around R ~ r. 

Examining the second integral in eq. ( |B-7| ), we may 
reasonably assume that the correlation Ba (R, r) be- 
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